La ville de Paris rend publiquement accès à de nombreux jeux de données sur la circulation routière de la ville. La découverte du site web Paris Data (“Paris Data” visité le 26.12.2021) nous a motivé à utiliser cette ressource dans notre projet de prédiction. La possibilité de relier les différents capteurs dans une structure de graphe nous a intéressé. D’une part, cela motive des études du réseau entier : peut-on construire un algorithme global, qui apprend les interactions dans le réseau entier ? D’e l’autre côté’autre part, l’examination de rues individuelles mène à des questions sur la corrélation locale de la circulation. Un grand boulevard se comporte-t-il comme ses avenues voisines ? Enfin, dans tout cela, nous nous poserons la question à quel horizon temporel nous réussirons à prédire la présence de voitures dans les rues de Paris.
Avant de nous lancer, nous présentons d’abord le jeu de données en détail, sa richesse mais aussi les difficultés qu’il apporte : Sa taille et les valeurs manquantes. Ayant passé une grande partie du projet sur cet aspect, nous présentons nos stratégies pour résumer et compléter ce data set.
Cela fait, nous formulons plusieurs problèmes de prédiction, à différents horizons temporels et en utilisons différentes parties du data set. (TODO: spécifier à la fin pour avoir une petite table de matières ici)
Sur le site Paris Data, on trouve énormément de jeux de données issus des activités de la ville, notamment sur le comptage routier (“Données Du Comptage Routier” visité le 26.12.2021). Ces données sont collectées en continu grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. Ces boucles et bornes de comptage mesurent principalement deux choses: le nombre de voitures qui passent et l’occupation du tronçon de route. Disponibles de 2014 jusqu’à aujourd’hui, ces deux données sont mesurées au pas horaire, résultant en environ \[ 8 \text{ ans } \times 365 \text{ jours } \times 24 \text{ heures } \times 3000 \text{ capteurs } = 2.1\times10^8 \text{ observations dans le data set.} \]
Chacune de ces lignes contient les informations suivantes, qui serons nos variables explicatives \(X_i\). Les noms des variables que nous utilisons dans notre analyse sont écrits en gras tout au long de cet exposé.
Les données historiques de circulation sont enregistrées sur le site web Paris Data sous le format .txt. Chaque fichier couvre toutes les bornes pour une semaine. Pour commencer, il faut donc télécharger les fichiers de toutes les semaines et les convertir en dataframes. Ce faisant, nous effectuons également un tri par identifiant de borne iu_ac.
Au vue de la quantité énorme de 200 millions d’observations, nous sommes obligés de réduire la taille du dataset. La première démarche consiste à garder les variables intéressantes: t_1h, nbCar, rateCar et iu_ac.
Ayant réduit la dimension de chaque observation, la deuxième étape est la réduction du nombre de capteurs de comptage. Pour ce faire, nous utilisons un modèle simplifié des rues de Paris. Etant donné que les capteurs couvrent les grandes axes routiers, dont le périphérique, notre objectif était d’agréger les données selon ces grands axes. Afin d’identifier où passent le plus de voitures, nous calculons les moyennes de nbCar pour chaque libellé, c’est-à-dire les moyennes à travers le temps et pour toutes les bornes associées à chaque libellé. En prenant les 200 libellés dont le nombre de voitures est le plus grand (hors périphérique), nous obtenons le graphique 2.1 qui représente très bien les axes routiers auxquels on s’attendait.
Ensuite, nous ignorons les libellés et regardons exclusivement ce graphique. De manière arbitraire, nous en déduisons le graphe de la figure 2.2, un schéma très simplifié de la circulation à Paris. Grâce à notre démarche d’utiliser les rues les plus fréquentées, cette abstraction nous permet d’avoir un modèle représentatif de la circulation parisienne.
Figure 2.1: Représentation en rouge des bornes d’observation associés aux 200 libellés les plus fréquentés (hors périphérique)
Figure 2.2: Graphe simplifié. Les points noirs indiquent les intersections entre nos arêtes.
Chaque arête dans le graphe 2.2 regroupe plusieurs capteurs dont nous agrégeons les données : après une collecte minutieuse et laborieuse des identifiants de toutes les bornes sur chacune des arêtes, nous associons pour chaque heure la moyenne de nbCar et rateCar des bornes de l’arête dont ils font partie. Cela a trois avantages:
Cette réduction des données résulte en la construction des jeux de données que l’on va manipuler et en la concrétisation du problème que l’on souhaite résoudre. Elle permet également d’obtenir des objets que nos ordinateurs sont capables de stocker efficacement. Il reste pourtant un problème: la présence de données manquantes.
En examinant des figures basiques des taux d’occupation, nous observons régulièrement une image comme dans le graphique 2.3 Malgré l’agrégation le long des arêtes, il reste un nombre considérable de valeurs manquantes, NA, dans le dataset. Afin d’exécuter des tâches de prédiction, il nous faut cependant des données complètes. Une partie conséquente de notre travail a donc été de compléter les valeurs manquantes.
Figure 2.3: Taux d’occupation pour une journée ouvrière avec une valeur manquante
Figure 2.4: Visualisation de valeurs manquantes pour 30 arêtes (gris = NA)
En regardant le graphique 2.4, on considère qu’il y a deux types de “trous” dans le data set:
Le fait que les arêtes ont de différentes structures de valeurs manquantes peut aussi être observé en regardant la répartition du pourcentage de NA’s à travers les dataframes dans la figure 2.5.
Figure 2.5: Distribution des valeurs manquantes
Nous observons que les données de taux d’occupation sont généralement plus complètes : ceci est dû au fait qu’il y a plus de capteurs de ce type. Par ailleurs, nous avons vérifié qu’une valeur manquante de nbCar implique toujours un NA dans rateCar.
Pour remplir les trous, la première étape consiste à ajouter des variables explicatives supplémentaires. En fonction du type de trou, nous verrons qu’ils auront un impact important sur la complétion.
Pour l’instant, nous disposons des variables nbCar et rateCar pour chaque arête et chaque heure entre 2014 et 2020. Nous exploitons deux propriétés structurelles des données pour obtenir de l’information supplémentaire:
D’abord, il s’agit de séries temporelles, donc nous rajoutons les valeurs historiques de nbCar et rateCar décalées d’une heure, d’un jour et d’une semaine (nbCarLaggedHour, nbCarLaggedDay, et nbCarLaggedWeek) comme variables explicatives. L’espoir étant que la circulation reste identique à travers ces cycles. Avec l’objectif spécifique de compléter les données manquantes, il est aussi pertinent d’ajouter ces valeurs du futur (rateCarFuturHour etc.), étant donné que le modèle d’imputation a pour but d’interpoler au lieu d’extrapoler. On fait l’hypothèse que toutes ces variables décalées seront efficaces pour remplir les trous locaux : une heure ou une journée manquante devrait facilement être inférée en utilisant des données qui sont proches dans le temps.
Pour les valeurs manquantes au milieu des grands trous, nous n’avons pas accès aux heures et jours d’avant car elles sont aussi manquantes. Cependant, grâce à l’interconnexion de nos données à travers le graphe présenté plus haut, les mesures de circulation d’une arête peuvent bien être expliquées en fonction de celle de ces voisines. Dans un effort manuel, nous avons pour chaque arête \(A\) collectionné des voisins \(V_i\), dans le sens que les voitures dans \(A\) passent ensuite par l’un des \(V_i\) ou bien l’inverse. Nous faisons ici encore des choix arbitraires de qui est voisin de qui pour deux raisons. Premièrement, inférer statistiquement quelle arête influence quelle autre nécessiterait des données déjà complètes pour une régression par exemple. Deuxièmement, utiliser toutes les autres 68 arêtes au lieu de juste 2 à 6 voisins augmenterait trop le temps d’exécution de la complétion.
Donnons un exemple illustrant le choix de voisins ainsi que la dénomination des nouvelles variables. Parmi les voitures qui entrent dans Paris par l’arête “pont amont - pont austerlitz,” il y a une partie considérable qui continue tout droit sur l’arête “pont austerlitz - chatelet.” Par conséquent, nous ajoutons la variable rateCar_pontausterlitzTOchatelet au dataframe de “pont amont - pont austerlitz.”
Jusqu’ici, nous n’avons utilisé que des informations déjà présentes dans le jeu de données. Dans la suite, nous ajoutons des variables extérieures qui pourraient également expliquer la circulation routière.
Variables temporelles
Le traffic routier étant relié à l’activité humaine, nous avons ajouté de nombreuses variables temporelles (notamment grâce à la librairie lubridate).
Index de la situation sanitaire en rapport avec le Covid-19
Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021). Cette variable permettra éventuellement une étude de l’année 2020 qui voit une circulation perturbé. Pour cela, il faudra utiliser des modèles qui s’adaptent très vite à l’influence de nouvelles variables.
Météo
A partir de données de l’Organisation Météorologique Mondiale (“Observations Météorologiques Historiques En France” visité le 26.12.2021), nous avons extrait 2 variables météorologiques : temperature la température en Kelvin et precipitation les précipitations dans les 3 dernières heures en mm. Ces données ont été mesurées à Athis-Mons en Essonne.
On dispose d’un relevé tous les 3 heures environ donc on procède à une interpolation pour compléter les données. Etant donné que les mesures sont uniformément réparties, on utilise une interpolation linéaire basique à l’aide de la fonction na.approx de la librairie zoo.
Avec un dataframe enrichi à notre disposition pour chacune des arêtes du graphe de circulation, nous pouvons passer au remplissage des valeurs manquantes. Nous sommes pourtant restreint dans notre recherche d’un algorithme d’imputation car plusieurs variables explicatives contiennent également des NA. D’une part, nous voulons remplir deux variables en même temps (nbCar et rateCar). D’autre part, les voisins n’ont pas moins de valeurs manquantes, laissant de grands trous dans le dataframe.
Heureusement, il y a un modèle qui peut être entraîné en dépit de valeurs manquantes parmi les variables : les forêts aléatoires constituées d’arbres CART. Une règle de décision dans un tel arbre peut être ignorée si la valeur de la variable nécessaire n’est pas renseignée. Le package miceRanger profite de ce fait en utilisant un algorithme dit d’imputation multiple : en commençant par la variable \(V_1\), une forêt aléatoire est entraînée sur les lignes où \(V_1\) n’est pas manquant. Puis, les autres lignes sont “prédites” par la forêt. Avec le nouveau dataset moins vide, le processus est reproduit et ainsi de suite (Van Buuren 2018). Nous nous arêtons pourtant à la deuxième itération vu qu’il n’y a pas d’intérêt à compléter les données des voisins.
Comme nous répétons le processus de complétion 69 fois, le temps d’exécution est de quelques heures sur nos ordinateurs. Nous n’avons donc pas le luxe d’optimiser des hyperparamètres comme nous ferons plus tard pour la tâche de prédiction. Le choix de 100 arbres par forêt et 7 variables considérées à chaque split semble un bon compromis. (voir le fichier preprocessing/missing_data.R sur git pour le code exact).
Au délà des capacités de complétion, les forêts aléatoires permettent de calculer des scores d’importance pour chaque variable explicatives. Ces scores peuvent aider à identifier lesquelles des variables ont le plus contribué à l’estimation de la cible (le score est calculé comme la réduction de variance à chaque split). Pour une arête du graphe, jussieu - saint michel, nous illustrons les variables qui ont le plus contribué à l’imputation :
Figure 2.6: Variables avec le score d’importance le plus haut pour rateCar de l’arête ‘’jussieu - saint michel’’
Figure 2.7: Variables avec le score d’importance le plus haut pour nbCar de l’arête ‘’jussieu - saint michel’’
Sur l’arête dont traitent les deux plots, on observe que rateCar et nbCar exhibent deux comportement assez différents. Pour le taux d’occupation, ce sont presque exclusivement les valeurs temporellement décalées qui contribuent à l’imputaion. Comme il y seulement 0.8% de valeurs à compléter, ceci confirme notre hypothèse que les petits trous d’une ou plusieurs heures sont très bien complétés par les valeurs des heures juste avant ou après. Pour le nombre de voitures, où la proportion de NA était plus importante, les deux variables ayant le plus réduit les variances aux splits sont des variables nbCar issues d’arêtes voisines. Les trous importants en nbCar entre Jussieu et Saint Michel semblent être le mieux expliqué par ce qui se passe autour. Cela confirme notre choix de variables pour combler les 2 types de trous. Un phénomène similaire peut s’observer dans les autres arêtes.(TODO Guillaume: Quoi d’autre écrire sur ceci ? Y a des scores intéressants dans RF (genre error Out of Bag) que j’aurais dû retenir de l’imputation ?)
Après l’imputation arête par arête, nous pouvons enfin compléter toutes les données grâce au fait que les variables reliées aux voisins ont été imputées dans les dataframes qui correspondent à ces voisins. Il est important de noter ici que nous allons dans la suite couper le jeu de données en deux (train - test) et que cette imputation sera effectuée individuellement sur chacune des deux parties. Faire ceci deux fois séparément est nécessaire pour préserver l’indépendance et permet ensuite d’estimer les erreurs de nos modèles.
Dans la suite, nous supposons que les données sont complètes. Il faut pourtant remarquer que ces forêts aléatoires font indirectement partie des modèles de prédiction que nous allons utiliser. De plus, nous gardons à l’esprit que nous avons introduit un biais dans nos données d’apprentissage comme de test.
Avec un jeu de données complet, nous pouvons enfin utiliser des algorithmes de Machine Learning pour faire de la prévision. Vu la complexité du jeu de données, il n’est même pas évident quelles questions se poser. Nous allons alors commencer par la définition de différents problèmes et de métriques de comparaison avant d’entraîner des modèles.
Il y a deux dimensions principales du problème que nous allons examiner dans la suite, en utilisant de différents modèles d’apprentissage:
Avant de nous lancer, nous précisons comment nous évaluons nos modèles et quels sont les benchmarks, les modèles naïfs, que nous voulons battre.
Train - test - split
On découpe notre jeu de données en 2 parties de proportion 2/3 et 1/3 : une partie apprentissage de 2014 à 2017 et une partie test de 2018 à 2019. Faire le découpage de cette manière assure qu’on est évalué sur les prédictions du futur plutôt que du passé. L’année 2020 est omise à cause du Covid, mais elle pourrait faire partie d’un autre projet sur le jeu de données en relation avec la pandémie !
Mesurer la performance
Pour évaluer les performances des modèles, nous considérerons le Root Mean Square Error (RMSE) qui représente la racine carrée du deuxième moment d’échantillonnage des différences entre les valeurs prédites et les valeurs observées. Une autre métrique souvent utilisée est le Mean Absolute Percentage Error (MAPE) qui n’est pas applicable à notre situation car un nombre important des valeurs de rateCar sont nulles et on ne peut pas diviser par elles.
Modèles naïfs à battre
Afin de pouvoir comparer nos modèles à une référence, une sorte de benchmark, on construit 3 modèles témoins. Ces modèles sont dit naïfs car extrêmement simpliste :
naiveModel1 prévoit le taux d’occupation d’une heure \(t \in \{ 0, \dots, 23 \}\) d’un jour \(j \in \{ 1, \dots, 7\}\) d’un mois \(m \in \{ 1, \dots, 12 \}\) d’une année \(a \in \{ 2018, 2019 \}\) en moyennant les taux d’occupation de l’heure \(t\) du jour \(j\) du mois \(m\) pour toutes les années \(a \in \{ 2014, \dots, 2017 \}\).
naiveModel2 fait de même en calculant la médiane et non la moyenne.
naiveModel3 est simplement le taux d’occupation à l’heure précédente rateCar_LaggedHour.
Pour ces trois modèles simples, on calcule leurs erreurs RMSE pour chacune des 69 arêtes pour les moyenner après (voir le fichier naive_model.r pour le code). Ceci nous donne :
| naiveModel1 | naiveModel2 | naiveModel3 |
|---|---|---|
| 5.844 | 5.899 | 3.955 |
Les 2 premiers modèles naïfs ont des scores similaires. Le troisième modèle naïf, qui consiste à prendre la variable rateCar de manière décalée, est meilleur. On interprète cette bonne performance par la très forte corrélation du taux de circulation d’une heure avec la suivante. Dans la suite, il sera intéressant de comparer ce modèle naïf avec des modèles plus complexes pour justifier leur utilisation.
Validation croisée
Pour la séléction de modèles et plus particulièrement pour la recherche de paramètres (d’un arbre CART par exemple), nous nous servirons beaucoup de la validation croisée. Ceci n’est pas évident pour les séries temporelles comme les données ne sont pas indépendant d’un ligne à l’autre. En plus, il vaut mieux éviter l’utilisation de données futurs pour prédire le passé. Entre ces problèmes et l’efficacité algorithmique, nous utilisons trois approches différentes:
Validation croisée sur blocs regroupant des mois successifs: on divise les 4 années de données d’apprentissage en blocs de 3 mois ou 6 mois (correspondant à 16 et 8 “CV-folds” respectivement), ce qui englobe les différentes saisonnalités de nos données. C’est la méthode qui nous semble le plus pertinent pour les séries temporelles car les blocs sont à peu près indépendants entre eux et on respecte les saisonnalités dans le découpage. Pourtant, il n’y a pas de package efficace dans R pour faire ceci de manière efficace et le fait de devoir le programmer manuellement chaque fois nous empêche de l’utiliser souvent.
Validation croisée progressive: on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante). Pour utiliser cette méthode, on utilise le package caret avec le paramètre timeSlice.
Figure 3.1: Validation croisée progressive (“Progressive Cross-Validation” visité le 22.02.2022)
Si nous ignorons la première idée numériquement plus dure, il semble plus naturel d’utiliser l’approche progressive pour notre type de données. Il y a pourtant le souci de largement surévaluer les données des premières deux années dans le score de CV, vu qu’elles sont représentées dans chaque couche. En plus, les différentes couches (folds) seront alors fortement corrélées comme elles ont plus de la moitie de leur contenu en commun. Les deux idées ont alors leurs failles.
Après tout, lors de nos recherches de paramètres, nous avons observé que les deux méthodes sélectionnent des paramètres assez similaires. La plupart du temps, nous nous contentons alors de simplement utiliser la validation croisée habituelle avec 16 couches. Ce chiffre est inspiré par l’approche de blocs sur mois successifs. L’avantage de ceci est que le package caret permet une implémentation efficace. Nous notons ici que nous fixons toujours une seed du générateur aléatoire pour le mélange des lignes effectué par caret. Ainsi le lecteur pourra relancer le code lui-même.
Commençons par la tâche qui est censée être la plus simple : Prédire la circulation dans les prochaines 60 minutes.
Au début, nous nous restreignons à une arête isolée dans le sens que nous ne prenons pas en compte la circulation sur ces voisins. Notre objectif initial était de faire une simple régression linéaire, mais il y a quelques relations pas linéaires du tout. Si on essaie par exemple, d’écrire le taux d’occupation d’une avenue en fonction de l’heure, le graphique suivant montre bien que les heures de pointe seraient mieux décrites par un polynôme que pour une droite :
Figure 3.2: Une journée typique sur l’arête
Pas conséquent nous travaillons plutôt sur des GAM que sur de simples régression. Constituer un modèle additif et en sélectionner le meilleur est une tâche qui peut être longue. C’est pourquoi nous adapterons une approche greedy pour la sélection de modèles où nous rajoutons une variable explicative à la fois. Idéalement ce sera celle qui apporte le plus de puissance prédictive supplémentaire (forward search).
Pour prédire une heure dans le futur, il est incontestable que l’heure juste avant contient le plus d’informations. La variable rateCar_LaggedHour sera donc la première à considérer. L’effet de l’occupation d’une heure à l’autre n’est pourtant pas homogène à travers les heures de la journée. En plus, comme nous voyons dans le plot supérieur gauche, il y a peut-être une relation linéaire entre les variables, qui a pourtant une grande variance, voir l’image à gauche en haut :
Les autres trois plots que nous venons de voir représentent la même comparaison rateCar_LaggedHours vs. rateCar, mais sur différentes heures. On a bien l’impression que des régressions linéaires pourront nous servir, vu que la variance qui est visiblement réduite. Comparer le nuage de points à 8 heures avec ceux de 12 et de 13 heures nous donne envie de regrouper certaines heures. Ce phénomène était déjà observable dans le graphique 3.2 : dans l’heure de pointe matinale (morning rushhour), la circulation augmente beaucoup entre 7 et 8 heures tandis qu’entre 12 et 13 heures elle reste plutôt constante.
Afin de regrouper les heures de la journées en périodes semblables, nous procédons de manière heuristique : pour chaque heure de la journée, nous effectuons une régression linéaire et comparons les coefficients (intercept et slope). Ceci donne le tableau suivant.
| hour | 0.000 | 1.000 | 2.000 | 3.000 | 4.000 | 5.000 | 6.000 | 7.000 | 8.000 | 9.000 | 10.000 | 11.000 | 12.000 | 13.000 | 14.000 | 15.000 | 16.000 | 17.000 | 18.000 | 19.000 | 20.000 | 21.000 | 22.000 | 23.000 |
| intercept | -0.148 | 0.038 | -0.218 | 0.039 | 0.007 | 0.163 | 0.414 | 1.155 | -0.630 | -1.185 | 1.691 | 2.253 | 1.688 | 2.278 | 0.756 | 1.778 | 2.442 | 1.527 | 1.230 | 1.529 | 2.020 | 0.709 | 2.252 | 1.675 |
| slope | 1.050 | 0.749 | 0.736 | 0.651 | 0.645 | 0.578 | 0.575 | 0.568 | 2.993 | 2.298 | 0.846 | 0.888 | 0.949 | 0.879 | 0.808 | 0.842 | 1.046 | 1.013 | 0.916 | 0.810 | 0.795 | 0.685 | 0.274 | 0.438 |
Les résultats sont congruents avec le graphique 3.2, par exemple l’explosion de slope à 8 heures où commence l’heure de pointe matinale. En nous basant sur ce tableau, nous regroupons de manière heuristique les heures en différents groupes: “night” de 0 à 6 heures, “7heures” à huit heures, “8heures” et “9heures” pareil, “noon” de 10 à 15 heures, “afternoon” de 16 à 17 heures, “evening_rush” de 18 à 20 heures et “evening” de 21 à 23 heures. Ainsi, chaque ligne du dataframe reçoit le nouveau attribut hour_groups. Avec ces groupes, nous pouvons établir notre premier modèle mod1, une régression linéaire avec un “mixed effect” hour_groups * rateCar_LaggedHour :
cluster_hours = function(h){
if (h %in% c(0,1,2,3,4,5,6,7)){return("night")}
if (h %in% c(7)){return("7heures")}
if (h %in% c(8)){return("8heures")}
if (h %in% c(9)){return("9heures")}
if (h %in% c(10,11,12,13,14,15)){return("noon")}
if (h %in% c(16,17)){return("afternoon")}
if (h %in% c(18,19,20)){return("evening_rush")}
if (h %in% c(21,22,23)){return("evening")}
else{return("asfd")}
}
df_train$hour_groups = as.factor(sapply(as.factor(df_train$hour), cluster_hours))
df_test$hour_groups = as.factor(sapply(as.factor(df_test$hour), cluster_hours))
# for (group in unique(df_train$hour_groups)){
# indices = df_train$hour_groups == group
# plot(df_train[indices,]$rateCar_LaggedHour, df_train[indices,]$rateCar,
# main = group, xlab = "rateCar_laggedHour", ylab = "rateCar")
# }
form1 = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour
mod1 = lm(data = df_train, formula = form1)
summary(mod1)
##
## Call:
## lm(formula = form1, data = df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.480 -0.836 -0.301 0.749 51.016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.846441 0.024047 35.20 <2e-16
## hour_groups8heures:rateCar_LaggedHour 2.201884 0.038810 56.73 <2e-16
## hour_groups9heures:rateCar_LaggedHour 1.935035 0.013804 140.18 <2e-16
## hour_groupsafternoon:rateCar_LaggedHour 1.097172 0.003656 300.10 <2e-16
## hour_groupsevening:rateCar_LaggedHour 0.572480 0.003810 150.26 <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.881992 0.002559 344.71 <2e-16
## hour_groupsnight:rateCar_LaggedHour 0.563815 0.011286 49.96 <2e-16
## hour_groupsnoon:rateCar_LaggedHour 0.928103 0.002810 330.23 <2e-16
##
## (Intercept) ***
## hour_groups8heures:rateCar_LaggedHour ***
## hour_groups9heures:rateCar_LaggedHour ***
## hour_groupsafternoon:rateCar_LaggedHour ***
## hour_groupsevening:rateCar_LaggedHour ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour ***
## hour_groupsnoon:rateCar_LaggedHour ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.421 on 35052 degrees of freedom
## Multiple R-squared: 0.8753, Adjusted R-squared: 0.8753
## F-statistic: 3.516e+04 on 7 and 35052 DF, p-value: < 2.2e-16
Toutes les variables ont une contribution significative d’après les tests de student. Calculons le score de validation croisée pour ce modèle simple pour obtenir
CV_gam = function(formula, df_train, no_folds){
K = no_folds
data = df_train
folds <- cut(seq(1,nrow(data)),breaks=K,labels=FALSE)
fold_scores = c()
start_time = Sys.time()
for(i in 1:K){
testIndexes <- which(folds==i,arr.ind=TRUE)
testData <- data[testIndexes, ]
trainData <- data[-testIndexes, ]
g = gam(data=trainData, formula=formula)
Y_test = testData$rateCar
Y_predict = predict(g, newdata=subset(testData, select = -c(rateCar)))
fold_scores = c(fold_scores, rmse(Y_test, Y_predict))
}
exec_time = Sys.time()-start_time
score = mean(fold_scores)
ret = data.frame(score = score, no_folds = no_folds, exec_time = exec_time,
formula = paste(as.character(formula)[c(2,1,3)], collapse = " "))
return(ret)
}
knitr::kable(CV_gam(form1, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.42125 | 8 | 1.119009 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour |
Afin de choisir la prochaine variable à rajouter pour expliquer rateCar, nous les explorons une par une :
form1_hour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, bs="cc", k = 24, by=weekendsIndicator)
mod1_hour = gam(form1_hour, data=df_train)
summary(mod1_hour)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour +
## s(hour, bs = "cc", k = 24, by = weekendsIndicator)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.101257 0.027548 39.98 < 2e-16
## hour_groups7heures:rateCar_LaggedHour 0.711020 0.135167 5.26 1.45e-07
## hour_groups8heures:rateCar_LaggedHour 2.304156 0.107999 21.34 < 2e-16
## hour_groups9heures:rateCar_LaggedHour 2.343870 0.033249 70.49 < 2e-16
## hour_groupsafternoon:rateCar_LaggedHour 1.026931 0.007463 137.61 < 2e-16
## hour_groupsevening:rateCar_LaggedHour 0.525550 0.006976 75.33 < 2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.819606 0.004874 168.18 < 2e-16
## hour_groupsnight:rateCar_LaggedHour 0.655250 0.016846 38.90 < 2e-16
## hour_groupsnoon:rateCar_LaggedHour 0.867199 0.004684 185.13 < 2e-16
##
## (Intercept) ***
## hour_groups7heures:rateCar_LaggedHour ***
## hour_groups8heures:rateCar_LaggedHour ***
## hour_groups9heures:rateCar_LaggedHour ***
## hour_groupsafternoon:rateCar_LaggedHour ***
## hour_groupsevening:rateCar_LaggedHour ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour ***
## hour_groupsnoon:rateCar_LaggedHour ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(hour):weekendsIndicator0 21.94 22 203.63 <2e-16 ***
## s(hour):weekendsIndicator1 21.44 22 75.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.894 Deviance explained = 89.4%
## GCV = 5.0041 Scale est. = 4.9966 n = 35060
knitr::kable(CV_gam(form1_hour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.2425 | 8 | 5.049379 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, bs = “cc,” k = 24, by = weekendsIndicator) |
#test_score = rmse(predict(mod1_hour, newdata = df_test), df_test$rateCar)
#print(test_score)
form1_weekdays = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(weekdays, bs="cc", k = 4, by=as.factor(winterHolidaysIndicator))
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_weekdays = gam(form1_weekdays, data=df_train)
summary(mod1_weekdays)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour +
## s(weekdays, bs = "cc", k = 4, by = as.factor(winterHolidaysIndicator))
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.859531 0.025096 34.25 <2e-16
## hour_groups7heures:rateCar_LaggedHour 0.853540 0.074459 11.46 <2e-16
## hour_groups8heures:rateCar_LaggedHour 2.186405 0.039039 56.01 <2e-16
## hour_groups9heures:rateCar_LaggedHour 1.928966 0.013898 138.80 <2e-16
## hour_groupsafternoon:rateCar_LaggedHour 1.095019 0.003699 296.01 <2e-16
## hour_groupsevening:rateCar_LaggedHour 0.569527 0.003870 147.16 <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.880178 0.002603 338.14 <2e-16
## hour_groupsnight:rateCar_LaggedHour 0.555410 0.011345 48.96 <2e-16
## hour_groupsnoon:rateCar_LaggedHour 0.925539 0.002894 319.86 <2e-16
##
## (Intercept) ***
## hour_groups7heures:rateCar_LaggedHour ***
## hour_groups8heures:rateCar_LaggedHour ***
## hour_groups9heures:rateCar_LaggedHour ***
## hour_groupsafternoon:rateCar_LaggedHour ***
## hour_groupsevening:rateCar_LaggedHour ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour ***
## hour_groupsnoon:rateCar_LaggedHour ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(weekdays):as.factor(winterHolidaysIndicator)0 1.948e+00 2 46.78 <2e-16
## s(weekdays):as.factor(winterHolidaysIndicator)1 5.889e-09 2 0.00 0.532
##
## s(weekdays):as.factor(winterHolidaysIndicator)0 ***
## s(weekdays):as.factor(winterHolidaysIndicator)1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.876 Deviance explained = 87.6%
## GCV = 5.8457 Scale est. = 5.8439 n = 35060
#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_weekdays, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.415 | 8 | 1.290045 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(weekdays, bs = “cc,” k = 4, by = as.factor(winterHolidaysIndicator)) |
form1_state = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(state, bs="cr", k = 5)
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_state = gam(form1_state, data=df_train)
#summary(mod1_state)
#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_state, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.42125 | 8 | 0.750025 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(state, bs = “cr,” k = 5) |
form1_toy = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(toy, bs="cr", k = 6)
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_toy = gam(form1_toy, data=df_train)
summary(mod1_toy)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour +
## s(toy, bs = "cr", k = 6)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.831406 0.024825 33.49 <2e-16
## hour_groups7heures:rateCar_LaggedHour 0.882487 0.074472 11.85 <2e-16
## hour_groups8heures:rateCar_LaggedHour 2.209279 0.038941 56.73 <2e-16
## hour_groups9heures:rateCar_LaggedHour 1.937016 0.013848 139.88 <2e-16
## hour_groupsafternoon:rateCar_LaggedHour 1.097766 0.003682 298.14 <2e-16
## hour_groupsevening:rateCar_LaggedHour 0.573190 0.003842 149.18 <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.882456 0.002584 341.51 <2e-16
## hour_groupsnight:rateCar_LaggedHour 0.563128 0.011324 49.73 <2e-16
## hour_groupsnoon:rateCar_LaggedHour 0.928841 0.002852 325.67 <2e-16
##
## (Intercept) ***
## hour_groups7heures:rateCar_LaggedHour ***
## hour_groups8heures:rateCar_LaggedHour ***
## hour_groups9heures:rateCar_LaggedHour ***
## hour_groupsafternoon:rateCar_LaggedHour ***
## hour_groupsevening:rateCar_LaggedHour ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour ***
## hour_groupsnoon:rateCar_LaggedHour ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(toy) 4.676 4.948 4.658 0.000219 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.875 Deviance explained = 87.5%
## GCV = 5.8582 Scale est. = 5.8559 n = 35060
#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_toy, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.42125 | 8 | 0.7370381 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(toy, bs = “cr,” k = 6) |
for (h in 0:23){
g = ggplot(data=df_train[(df_train$hour == h),], aes(y = rateCar, x = nbCar_LaggedHour))+
geom_point() + xlab("nbCar_laggedHour") +
ylab("rateCar (présent)") +
labs(title = paste0("Heure ", h))
plot(g)
}
form1_nbCarHour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(nbCar_LaggedHour, bs="cr", k = 6, by=hour_groups)
mod1_nbCarHour = gam(form1_nbCarHour, data=df_train)
summary(mod1_nbCarHour)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour +
## s(nbCar_LaggedHour, bs = "cr", k = 6, by = hour_groups)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.845837 0.082954 22.251 < 2e-16
## hour_groups7heures:rateCar_LaggedHour 0.905683 0.202185 4.479 7.51e-06
## hour_groups8heures:rateCar_LaggedHour 1.560141 0.162439 9.604 < 2e-16
## hour_groups9heures:rateCar_LaggedHour 1.361042 0.066729 20.396 < 2e-16
## hour_groupsafternoon:rateCar_LaggedHour 0.929243 0.009734 95.468 < 2e-16
## hour_groupsevening:rateCar_LaggedHour 0.536947 0.010040 53.481 < 2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.818969 0.005512 148.568 < 2e-16
## hour_groupsnight:rateCar_LaggedHour 0.455403 0.025760 17.679 < 2e-16
## hour_groupsnoon:rateCar_LaggedHour 0.804931 0.006956 115.710 < 2e-16
##
## (Intercept) ***
## hour_groups7heures:rateCar_LaggedHour ***
## hour_groups8heures:rateCar_LaggedHour ***
## hour_groups9heures:rateCar_LaggedHour ***
## hour_groupsafternoon:rateCar_LaggedHour ***
## hour_groupsevening:rateCar_LaggedHour ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour ***
## hour_groupsnoon:rateCar_LaggedHour ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(nbCar_LaggedHour):hour_groups7heures 3.610 3.884 10.44 <2e-16 ***
## s(nbCar_LaggedHour):hour_groups8heures 4.083 4.374 95.04 <2e-16 ***
## s(nbCar_LaggedHour):hour_groups9heures 4.086 4.277 159.35 <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsafternoon 4.985 5.000 36.38 <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsevening 4.847 4.981 54.04 <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsevening_rush 4.761 4.949 17.89 <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsnight 3.250 3.682 132.29 <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsnoon 4.797 4.979 27.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.885 Deviance explained = 88.5%
## GCV = 5.4364 Scale est. = 5.4297 n = 35060
#plot(mod1_nbCarHour, residuals=T, rug=T, se=T, pch=20)
#draw(mod1_nbCarHour, residuals=T)
knitr::kable(CV_gam(form1_nbCarHour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")
| score | no_folds | exec_time | formula |
|---|---|---|---|
| 2.4525 | 8 | 3.326821 secs | rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(nbCar_LaggedHour, bs = “cr,” k = 6, by = hour_groups) |
Après cet exercise de bricolage avec GAM, nous passons à d’autres modèles dont le tuning pourra se faire de manière plus automatique. Pour commencer, nous nous intéressons à des arbres CART et nous nous servons du package rpart pour les implémenter. Au lieu d’examiner la contribution de chaque variable, le choix des \(k\) etc., il suffit de tuner quelques paramètres de l’algorithme. Vu que nous ferons un tuning étendu pour les forêts aléatoires après, nous nous focalisons ici sur un seul paramètre qui est de plus haute importance : le “complexity parameter” cp. Dans la construction de l’arbre, si le gain de performance après une découpe n’est pas meilleur d’un facteur de cp, alors la découpe n’est pas réalisée. Plus cp est grand, moins l’arbre sera alors complexe, dont le nom.
Pour optimiser cp, nous utilisons de la validation croisée et une recherche sur deux grilles : la première de 0.0 à 0.1 et la seconde, plus fine, de 0.0 à 0.01.
Nous affichons les scores de CV obtenues (voir le fichier 1_heure/tree.R pour le code) :
Figure 3.3: Recherche du paramètre cp optimal
La valeur optimale d’après cette CV est alors cp=0, c’est-à-dire que l’on accepte toutes les découpes et que notre arbre est un arbre profond. On maintient à leur valeur par défaut les autres paramètres de CART et on entraîne le modèle sur notre train dataset. Nous faisons ceci dans deux cas différents qui représentent deux dimensions locales différentes : rappelons-nous que pour la complétion, nous avions rajouté des informations sur les voisins pour enrichir les données. Pour voir si ces mêmes informations nous servent aussi maintenant, nous entraînons un arbre avec et un CART sans données voisines. Voici les scores sur le test set dans les deux cas :
| Arbre avec voisins | Arbre sans voisins |
|---|---|
| 3.783 | 3.582 |
Premièrement, on note dans les deux cas une performance légèrement meilleure que le modèle naïf qui recopie l’heure d’avant. Ce qui est pourtant frappant est que le rajout d’information sur les voisins baisse le score, indiquant qu’il y a soit de l’overfitting, soit un problème de grande dimension qui perturbe le bon fonctionnement des arbres aléatoires. En tout cas, ce résultat sème le doute sur l’utilité de la composante spatiale de nos données. Mais peut-être, il faut juste utiliser des techniques plus puissantes, comme les forêts aléatoires qui rassemblent plusieurs arbres en utilisant du bagging.
Pour cette méthode d’ensemble basé sur les arbres de Breimann, le package ranger nous servira à construire les forêts et le package vip aidera à visualiser l’importance des variables. Comme pour les arbres, il y a des paramètres à optimiser avec de la CV pour que la forêt apprenne bien mais pas trop la circulation parisienne. Nous nous focalisons sur deux paramètres que nous optimisons l’un après l’autre. C’est une approche moins prometteuse qu’une recherche sur grille, mais elle diminue le temps de calcul et nous permet de visualiser la recherche de paramètres :
min.node.size (minimal node size). L’algorithme ne découpe pas une feuille dont le nombre d’éléments est inférieur ou égal à sa valeur. Comme cp avant, ce paramètre contrôle la taille des arbres.Figure 3.4: Recherche du paramètre min.node.size optimal
On remarque que l’erreur décroît jusque la valeur 1 avec un pic négatif à 9. Choisir min.node.size=1 donnerait des arbres de profondeur maximale, rendant le modèle susceptible à sur-apprendre. Par conséquent nous nous contentons de choisir 9, malgré les fluctuations importantes des CV-scores dans le graphique à droite. Ces flucuations pourraient être dues au fait qu’ici on a fait une validation croisée habituelle (voir nos explications sur CV en haut) qui mélange les lignes corrélées avant de les séparer en blocs. Dans un cadre iid, nous nous attendrions à la forme typique d’un “coude” par exemple.
Figure 3.5: Recherche du paramètre mtry optimal
Celle-ci est l’unique fois où nous affichons la méthode de validation croisée par timeslice. On remarque que l’allure globale des courbes sont identiques mais que le coude est plus lisse avec timeslice. On choisit de poser \(\textbf{mtry} = 10\), qui semble être un bon choix d’après les deux graphiques.
Le dernier paramètre à déterminer est ntree, le nombre d’arbres dans la forêts. Ce n’est pas un paramètre à optimiser comme les 2 précédents car le sur-apprentissage dans les forêts aléatoires n’est pas lié au nombre d’arbres. On optimise ntree afin de déterminer à partir de quel nombre d’arbre on gaspille du temps d’exécution pour un faible gain de performance.
Figure 3.6: Recherche du paramètre ntree ‘’optimal’’
On voit que le temps d’exécution est linéaire selon le nombre d’arbres et que la courbe de l’erreur marque un coude à partir de \(150\), valeur que l’on choisit pour **ntree*.
Avec ces trois paramètres choisis, nous entraînons la forêt aléatoire et l’évaluons sur le test set. Comme avec les arbres, nous faisons ceci deux fois, une fois avec voisins et une fois sans.
| Forêt aléatoire avec voisins | Forêt aléatoire sans voisins |
|---|---|
| 3.092 | 3.129 |
Comme attendu, les performances des forêts aléatoires sont généralement meilleures que celles des arbres aléatoires. Pourtant, contrairement à ces derniers, les arêtes voisins semblent avoir un effet plutôt positif. Si nous nous intéressons davantage au rôle des voisins, nous pouvons exploiter le fait que les forêts aléatoires livrent des scores d’importance pour chaque variable. Sur l’arête pont amont - pont austerlitz, par exemple, le ranking d’importance des variables est le suivant :
Figure 3.7: Scores d’importance dans l’ordre décroissante (15 premières arêtes) pour l’arête ‘’pont amont - pont austerlitz’’
Sur cette arête, les voisins dominent même les scores d’importance. Ceci nous donne une motivation pour davantage explorer le rôle des voisins et ainsi le rôle de la dimension spatiale du problème de prédiction de la circulation :
Une question immédiate est quelles arêtes profitent le plus des données de leurs voisins. Nous essayons de visualiser ça sur le graphique suivant où pour chaque arête du graphe routier, nous visualisons la différence des erreurs test avec et sans voisins. Le trait vertical noir pointillé marque la séparation entre les rues à l’intérieur Paris et les parties du périphérique. Ces dernières se trouvent à droite de la ligne. Le trait horizontal rouge pointillé est un seuil arbitraire à partir duquel nous interprétons que l’ajout des voisins impacte la performance de la forêt significativement. Les arêtes qui dépassent le seuil en bas (et souffrent donc de l’inclusion de leurs voisins) sont représentées par une croix rouge, les autres par une croix noire.
TODO Jakob discuter avec Guillaume. Seules 12 forêts sur 69 ont une performance ‘’significativement’’ meilleure lorsqu’on ajoute les voisins. De plus, pour les arêtes du périphérique, on remarque que les différences sont plus importantes que pour Paris, qu’inclure les voisins est néfastes (à 2 exceptions près) et qu’en général, la circulation sur ces arêtes est moins bien prédite. En effet, les scores pour Paris sont 2.747 (avec voisins) et 2.898 (sans voisins), et les scores du périphérique sont 4.069 (avec voisins) et 3.787 (sans voisins).